#include <stdio.h>
#include <iostream>
#include <fstream>
#include <stdlib.h>
#include <string>
#include <gsl/gsl_math.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>

#include "IOParams.h"
#include "Gas.h"


using namespace std;
int main(int argc, char **argv){
	// TODO:Change initial conditions
  IOParams io(argc, argv);
  Gas g_old(io.N);
  Gas g_new(io.N);

  // Initiate angles and velocities.
	g_old.initCoord2();
	g_old.initVeloc2(io.T1);  
	

  ofstream ofile, ofile2;
	ofile.open((io.outfilename).c_str());
	ofile2.open((io.outfilename2).c_str());
	
  /**** ----------- Begin simulation   ----------- ****/
	// 1. Calculate ___force_old___.
	g_old.calculateVec_m();	  			//--  1.1 Calculate m_old.
	g_old.calculateVec_M();					//--  1.2 Calculate M_old.  
	g_old.calculateForce(io.eps);		//--  1.3 Calculate force_old with M_old

  while (io.ttime < io.total_time){
  	if (io.ttime % 2 == 0){		
			// 1. Calculate ___force_old___.
			// XXX: Was calculated in a previous step
			
			// 2. Calculate ___theta_new___
			g_new.updateCoord(&g_old, io.dt, io.dt2);
				 
			// 3. Calculate ___force_new___
			g_new.calculateVec_m();					//--  3.1 Calculate m_new
			g_new.calculateVec_M();					//--  3.2 Calculate M_new
			g_new.calculateForce(io.eps);		//--  3.3 Calculate force_new with M_new
				
			// 4. Calculate ___omega_new___
			g_new.updateVeloc(&g_old, io.dt);
			
			if (io.ttime % io.taux == 0){
				// Calculate system energy.
				g_old.calculateEkin();
				g_old.calculateEpot(io.eps);
				// Save data.	
				g_old.writeMacroState(&ofile, io.ttime, io.dt);
				g_old.writeMicroState(&ofile2, io.ttime, io.dt);
			}
		// -------------- End EVEN	---------------------------
  	}else{
  		// XXX: WARNING this is an INTERCHANGE
  		// 1. Calculate ___force_old___.
			// XXX: was calculated in the previous step
			
			// 2. Calculate ___theta_new___
			g_old.updateCoord(&g_new, io.dt, io.dt2);
				 
			// 3. Calculate ___force_new___
			g_old.calculateVec_m();					//--  3.1 Calculate m_new
			g_old.calculateVec_M();					//--  3.2 Calculate M_new
			g_old.calculateForce(io.eps);		//--  3.3 Calculate force_new with M_new
				
			// 4. Calculate ___omega_new___
			g_old.updateVeloc(&g_new, io.dt);
			
			if (io.ttime % io.taux == 0){
				// Calculate system energy.
				g_new.calculateEkin();
				g_new.calculateEpot(io.eps);
				// Save data.	
				g_new.writeMacroState(&ofile, io.ttime, io.dt);
				g_new.writeMicroState(&ofile2, io.ttime, io.dt);
			}
			
		// -------------- End ODD	---------------------------
  	}
  	io.ttime++;
  }
  cout<<endl;
		
  /****----------- End simulation   ----------- ****/  
	ofile.close();
  ofile2.close();
  
}

